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We analyze free expansion of a trapped one-dimensional Bose gas after a sudden release from the confining 
trap potential. By using the stationary phase and local density approximations, we show that the long-time 
asymptotic density profile and the momentum distribution of the gas are determined by the initial distribution 
of Bethe rapidities (quasimomenta) and hence can be obtained from the solutions to the Lieb-Liniger equations 
in the thermodynamic limit. For expansion from a harmonic trap, and in the limits of very weak and very strong 
interactions, we recover the self-similar scaling solutions known from the hydrodynamic approach. For all other 
power-law traps and arbitrary interaction strengths, the expansion is not self-similar and shows strong depen¬ 
dence of the density profile evolution on the trap anharmonicity. We also characterize dynamical fermionization 
of the expanding cloud in terms of correlation functions describing phase and density fluctuations. 
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The vast majority of natural and laboratory-induced phe¬ 
nomena occur in interacting many-particle systems that are 
away from the equilibrium. Yet the nonequilibrium dynamics 
of such systems remains an unsolved problem both in quan¬ 
tum physics and some areas of classical physics such as fluid 
dynamics. Examples include quantum and classical turbu¬ 
lence, dynamics across phase transitions, and plasma insta¬ 
bilities, to name a few. Ultracold atomic gases have recently 
emerged as a particularly promising platform to gain new in¬ 
sights into aspects of nonequilibrium dynamics of quantum 
many-body systems in El and, more generally, into nonequi¬ 
librium statistical mechanics. In particular, there has been a 
surge of research activity in the study of the dynamics af¬ 
ter a sudden quench of the system’s parameters, exploring 
the mechanisms of relaxation and the role of integrability in 
approaches to equilibrium EHn (for further references, see 

iniEi). 

In this Letter, we study far-from-equilibrium behavior of 
a trapped quantum gas after a sudden quench of the confin¬ 
ing potential. More specifically, we investigate free expan¬ 
sion of an interacting one-dimensional (ID) Bose gas instanta¬ 
neously released from the confining trap potential V(x). This 
is a paradigmatic example of a “quantum explosion” problem 
in an experimentally realizable system that can be described 
by an integrable microscopic model—the Lieb-Liniger model 
ca of delta-interacting bosons in one dimension. The exact 
integrability of the model offers an opportunity to investigate 
the expansion dynamics using theoretical methods that would 
have otherwise been inapplicable. At the same time, integra¬ 
bility implies that the underlying system lacks any mechanism 
of thermalization, which in turn poses a question of applica¬ 
bility of the standard hydrodynamic approach that was previ¬ 
ously used to describe the dynamics of this system lEHISl. 
Combining these aspects together gives us a unique opportu¬ 
nity to: (i) solve the quantum explosion problem in a nontriv¬ 
ial manner, and (ii) benchmark the predictions of the hydro- 
dynamic approach against those obtained here. 

In this work, we treat the simplest case of expansion from 


the zero-temperature ground state of the trapped gas. As we 
show below, the asymptotic density and momentum distribu¬ 
tions of the gas after a sufficiently long expansion time (once 
the expansion becomes ballistic) can be obtained from the 
initial distribution of quasimomenta of the trapped (nonuni¬ 
form) gas using the stationary phase approximation. This pro¬ 
motes the initial quasimomentum distribution from an auxil¬ 
iary quantity—which has so far only been used to derive ther¬ 
modynamic quantities—to the status of an observable phys¬ 
ical property. The initial quasimomentum distribution itself 
is calculated by combining the exact solutions of the uniform 
Lieb-Liniger model and the local density approximation. 

To start, we consider Hamiltonian evolution in free space 
of the many-body wave function, which—immediately prior 
to the removal of the trap potential at t = 0—describes the 
ground state of N trapped particles and is expanded in terms 
of the eigenfunctions of the uniform Lieb-Liniger model. 


T'(a:i,. ..,XN;t) = 


(27r)^/2 


dki ... dkf4 




Here, b{ki,... ,k]s[) are the expansion coefficients of the 
initial wave function ll20l . which depend on N differ¬ 
ent quasimomenta {kj} (also referred to as rapidities, 
and having units of wave numbers) and are normalized 
to jdki...dkN \b{ki,... ,k]s[)\'^ = 1- The phase 

0(fci,..., kjsf) = J2j<i ta,n~^[h‘^{ki — kj)/mg\ arises from 
two-body collisions by the diffractionless delta-function inter¬ 
action potential with the strength g ED- 

The expansion coefficients b{ki,... ^kjy) determine the 
joint W-particle probability distribution of quasimomenta 
\b{ki ,..., fcAr)p. In the single-particle sector, i.e., after in¬ 
tegration over all quasimomenta but one, they give the quasi¬ 
momentum distribution of the trapped gas, 

g{k) = N J dk 2 ...dkN\b{k,k 2 ,...,kN)\‘^, (2) 

with the normalization jg{k)dk = N. 
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To proceed, we note that the only time-dependent term in 
the integrand of Eq. 0 is the phase of the last exponen¬ 
tial. This exponential will, for sufficiently long times, develop 
fast oscillations as functions of kj compared to the remain¬ 
ing time-independent terms. Therefore, in the long-time and 
large-distance limit, the asymptotic form of the wave function 
can be simplified significantly by using the stationary phase 
approximation ll22l : the main contribution to the integral in 
Eq. Q comes from the stationary phase points EOll satisfying 


hkj Xj 
m t 


(3) 


and leading to the following asymptotic wave function; 



{mxi 

mxN\ 

\ht) ' 

K ht ’ 

’’ ht ) 


xe 




(4) 


The corresponding asymptotic density distribution 
Pac,{x,t) = N JdX 2 ■ ■ ■dXN\'i’ac{x,X 2 , ■ ■ ■ ,XN;t)\'^ Can 
therefore be found, using Eq. Q, as 



with Jp^{x,t)dx = N. Thus, the density profile after long 
expansion time is determined by the rescaled shape of the ini¬ 
tial quasimomentum distribution g{k)', finding this distribu¬ 
tion constitutes, therefore, a key task of the present work. 

The asymptotic wave function 0 also determines the 
asymptotic momentum distribution of the Indeed, the 
Eourier transform 'koo(fcii ■ ■ ■ of Eq. W is again dom¬ 

inated by the stationary phase points satisfying Eq. Q, re¬ 
garded now as conditions on the positions Xj. The result is 

^^ik„ ..., kN;t) = {-ifb{ku ..., ( 6 ) 

Integrating |4'oo(A:i,..., fc^r; over all momenta but one 
and using Eq. (j^, one obtains the asymptotic momentum dis¬ 
tribution, 

noo{k,t) = g{k), (7) 





FIG. 1. Illustration of the evolving local quasimomentum distribution 
f{k, x;t) = f{k, x-hkt/m-, 0), where Vu^a.yL = hk^^^/m. 


Einding the coefficients b{ki, fc 2 • ■ ■, fcw) and hence the dis¬ 
tribution g{k), Eq. 0, is equivalent to solving for the ground 
state of the (nonintegrable) trapped gas and thus constitutes 
a formidable task for systems with large N ||23|. However, 
as we argue here, for large N, the quasimomentum distribu¬ 
tion g{k) can be found approximately within the local density 
approximation (EDA) ll24l . 

The EDA is invoked by assuming that the initial trapped 
cloud can be divided into subsystems of length AL small 
compared its overall size R, so that the density p{x) = p{x,t = 
0) within each subsystem centered at x is approximately con¬ 
stant. At the same time AL has to be sufficiently large com¬ 
pared to the microscopic correlation length so that the lo¬ 
cally uniform subsystem can be treated via the solution of 
the Lieb-Liniger integral equation HI in the thermodynamic 
limit. Detailed conditions for the applicability of the EDA 
to trapped ID Bose gases have been discussed in Ref. ll^ . 
in addition to being verified experimentally ll26l . Generally 
speaking, the EDA is expected to be very good in the bulk of 
the atomic cloud for sufficiently large R (much larger than the 
length scale associated with the trapping potential), breaking 
down only in the small vicinity (~ ^q) of the cloud edge. 

The solution to the Lieb-Liniger integral equation for each 
region gives the local quasimomentum distribution /(fc, x;t = 
0) ifTrl corresponding to density p{x), obtained for the local 
value of the chemical potential p{x) =pq — V (x) ll25l . where 
Po is the global chemical potential. Integrating f{k,x;t = 
0) over X will give the quasimomentum distribution of the 
trapped gas in the local density approximation. 


implying that the initial quasimomenta of the trapped gas are 
mapped to real momenta of the expanded cloud 19] [221. Then, 
the result of Eq. 0 for the density profile can simply be 
viewed as a consequence of the ballistic position-momentum 
correlations, Eq. established in the long-time asymptotic 
regime, after the interaction energy has converted into the ki¬ 
netic energy of expanding particles. 

The characteristic expansion time te ensuring the applica¬ 
bility of the stationary phase approximation can be estimated 
by requiring that the fastest particles, moving with the velocity 
Umax Co, where cq is the speed of sound, overtake all slower 
particles. This is equivalent to cgte becoming larger than the 
characteristic size R of the initial cloud, and hence t^^R/cQ. 
Eor expansion from a harmonic trap with frequency wq this 
yields fe ~ 1 /wo in the Thomas-Eermi approximation, in both 
the weakly and strongly interacting regimes (see below). 


g{k) = J f(k,x;0)dx, (8) 

whereas p{x) = f f{k, x; 0)dk, with the normalization condi¬ 
tion J f{k,x\0)dkdx = N. 

It is insightful to see how the asymptotic density evolu¬ 
tion Eq. 0 can be obtained directly from the local quasimo¬ 
mentum distribution f{k,x-,t) if the latter is provided with 
a semiclassical time dependence of the form f{k,x;t) = 
f{k,x — hkt/m;0), for t ^ tg. This choice makes use of 
the ballistic expansion relationship x = hkt/m of Eq. 0 and 
is illustrated in Eig. It is clear that such an evolution leaves 
the (quasi)momentum distribution intact, 

g{k, t)= f{k, X — Hkt/m; 0) dx = g{k), (9) 
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whereas the density distribution evolves according to 

p{x^ t) dk =Jf{k, X — Hkt/m; 0) dk. (10) 


Introducing a new variable y = x — hkt/m, Eq. (10 1 can be 
rewritten as p{x,t) = Jf dy, where we 

can further neglect my/ht in the first argument as the main 
contribution to the integral comes from the values of y that 
are of the order of the initial size of the cloud R and values 
of a; oc f that are much larger than R in the long-time limit. 
Using Eq. one then obtains the same result as in Eq. 0 , 
t) = ^J dyf {^,y-0) as anticipated. 

We now apply our approach to expansion from power-law 
traps, V{x) = ^a,^\x\'', where i/ > 2 and is the confine¬ 
ment strength (in the important case of a harmonic trap, i/ = 2 
and a 2 = mcug). The problem can be treated analytically in 
two limiting cases: a weakly interacting gas (70 <C 1 ) and a 
strongly interacting gas in the Tonks-Girardeau (TG) regime 
(70 —>■ 00 ), where 70 = p^) is the dimensionless in¬ 

teraction strength in the trap center, with po being the peak 
density of the initial trapped sample. The intermediate regime 
can be addressed by finding numerically the local quasimo¬ 
mentum distribution f{k,x;0) via the solution of the Lieb- 
Liniger integral equation IH and then using Eqs. (j^ and Q. 
The results of such a numerical treatment for a harmonic trap 
are shown in Eig. (a), whereas the analytic results (see be¬ 
low) for v^2 are illustrated in Eigs.j^b) and 2(c). 

In the weakly interacting regime (70 <C 1), the local semi- 
classical distribution of a trapped gas with the density p(x) is 
given by Cl 


tn 1 4:H^p{x) 

f{k,x;0) = - 

ZTT y mg 


MiL2 

for \k\<K{x), (11) 


and f{k,x;0) = 0 otherwise, with K(x) = y/Amgp{x)/. 

The equation of state for a uniform gas in this regime is p = 
p/g, and therefore the density profile of the trapped sample in 
the Thomas-Eermi limit is given by 

p(x) = pix)/g = po (1 - \x\‘'/R''), for |x|<i?, ( 12 ) 

and p{x) = 0 otherwise. Here, R = (2pqI 
is the Thomas-Eermi radius, and po = Po/i? = [(1 + 
vYa.vN'" is the peak density found from 

the normalization condition iV=Jp(x)(ia; = 2 i?poi^/(l + v)- 
Integrating the distribution function /(/c,x;0), Eq. (Ill, 
over position gives 

) 7 „3, 


TT \ mg J Amgpo J 


for |/c| < y/AmgpoJ^, and g{k) = 0 otherwise. Here, R = 

fo dy (1 - = (2-^^a/2+\/^.) ^ with r(z) being the 

gamma function. The asymptotic density distribution, from 
Eq. (0, is then determined by 


Poo (Xj f) — 


4R 


Po 


TT X{t) 


1 - 


2 \ 1/2-tl/i^ 


A(f)2i?2 




(14) 


(a) 


(b) 


(c) 




s 



-1 -0.5 0 0.5 1 

x/R, k/Q 


FIG. 2. (Color online) Examples of the initial density profile p(x) 
(black solid lines) and the quasimomentum distribution g{k) (blue 
dash-dotted lines), which determines the shape of the asymptotic 
density profile poo{xR) and the momentum distribution noo{k,t), 
Eqs. and with Q being the maximum (quasi)momentum. (a) 
Expansion from a harmonic trap, with the solid and dash-dotted lines 
corresponding to the numerical results for 70 = 1. The analytic re¬ 
sults (coinciding with the hydrodynamic self-similar solutions) in 
the weakly (70 <C 1 ) and strongly (70 —> 00 ) interacting regimes 
are shown, respectively, by the dashed-magenta and dotted-red lines. 
The numerical results in these regimes, obtained for 70 = 2.5 x 10“"^ 
and 70 = 200 , are indistinguishable from the respective analytic 
curves and are omitted from the graphs for clarity, (b) Main curves 
are for a highly anharmonic trap with v = 14 El and 70 ^ 1 ; the 
semicircle (dashed red line) corresponds to the limiting behavior of 
g{k) for a box potential (i/ —>■ 00 ). (c) Strongly interacting regime, 
for the same v = 14. 


where we have introduced a dimensionless parameter A(t) = 
t/te = 2cot/R, with Co = yjPol 1 x 1 being the sound veloc¬ 
ity in the trap center. By comparing this result with the initial 
density distribution, Eq. it is now easy to see that for 
u = 2, for which /2 = 7 r /4 and te = 1/ {y/2uJo), we immedi¬ 
ately reproduce the scaling solution of Refs. GSIIISI, in which 
A(f) = y/2ujQt takes the meaning of the single scaling param¬ 
eter. In the hydrodynamic approach, \{t) is obtained from the 
scaling equation A = Wq/A^ EOll . We can also immediately 
conclude that such a self-similar scaling solution is not sup¬ 
ported by any other power-law trap potential. In particular, 
the case of 1 / = 14 illustrated in Eig. 0(b) shows a dramatic 
difference between the density profiles of the initial and ex¬ 
panded clouds {cf. 1281 ). 
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In the TG regime (70 —00), the local semiclassical distri¬ 
bution f{k, x; 0) is given by l[T4l 

f{k,x;0) = l/27r, for |fc| < Trp{x), (15) 


and f{k,x;Q) = 0 otherwise, where 7rp{x) is the maximum 
quasimomentum coinciding with the Fermi momentum of an 
ideal uniform Fermi gas of density p{x). 

The density profile p{x) is found from the equation of state 
of a uniform system p = y^2mp ./yielding 


p{x) = y'2mp{x)/{hny = po\/l - \x\''/R'', (16) 


for |a;| <R, and p{x) = 0 otherwise. Here, R = (2p.Qj^ 
is the Thomas-Fermi radius and po = = 

is the peak density, with be¬ 
ing the same numerical coefficient as in Eq. 

Integrating f{k, x; 0), Eq. (151, over position gives 


g{k) = 3 ( 1 - 


7^ Vo 


Xjv 


(17) 


for I/cI < 7rpo, and p(fc) = 0 otherwise. The asymptotic den¬ 
sity distribution is therefore given by 


Poo{x,t) 


Po f. _ 

X{t) \ X{tyR^ 




( 18 ) 


where X{t) = t/te = cot/R and cq = ^/2poJm is the sound 
(Fermi) velocity in the trap center. By comparing pao{x,t) 
with the initial density distribution, Eq. ( [T6] l, we immedi¬ 
ately see that, for finite ly (see also 1291 ), a self-similar scal¬ 
ing solution is again supported only by a quadratic potential, 
in which case A(t) = uJot. In the hydrodynamic approach, 
this asymptotic behavior is obtained from the scaling equa¬ 
tion X = lUq/X^ EOll (see also iQOl ). which also follows from 
the exact treatment of Ref. ED. 

Considering now the coherence properties of an expand¬ 
ing ID Bose gas, we note that the only length scale entering 
into the asymptotic momentum distribution rioo {k, t) [through 
Eqs. Q and either ( [T3] l or ( [T7) l] and hence into the respec¬ 
tive one-body density matrix {x, x'] t) is the microscopic 
correlation length ^0 = fi/^/mpo i25][32l, corresponding to 
the healing length ^0 = ^1 y/mgpQ for the weakly interacting 
gas and the mean interparticle separation ^0 = 1/po for the 
TG gas. In all cases, ^0 is much smaller than the size of the 
sample R, which implies complete loss of phase coherence (if 
there was any initially) typical of fermions and can be viewed 
as a manifestation of “dynamical fermionization” discussed 
in Refs. EIlESl. Such a loss of phase (or first-order) coher¬ 
ence with expansion, which we note does not follow from the 
hydrodynamic approach, is indeed the case for a weakly in¬ 
teracting gas, for which the initial (zero-temperature equilib¬ 
rium) coherence length ll34l is exponentially 

large and can typically be much larger than R. 

Another manifestation of dynamical fermionization during 
expansion can be seen in the asymptotic behavior of the same- 
point two-body correlation function g^'^^x, x;t): it acquires 


EOll a scaling cx 1/t^, implying suppressed correlation in the 
long time limit. Such a suppression indicates dynamical ap¬ 
proach to the fermionized TG regime, where g^"^^ {x, x) = 0 
due to an effective Pauli exclusion Esnsa. Moreover, our 
dynamical result can be written as g^'^\x,x',t) oc l/^(x,t)^ 
using the inverse scaling of the instantaneous interaction con¬ 
stant ^{x,t) = mg/h?pao{x,t) with density pao{x,t) oc 1/t. 
Such a scaling of the -function with 7 is indeed typical of 
an equilibrium TG gas E2I361. It must be noted though that 
this result is by no means an indication of equilibration dur¬ 
ing expansion as the time scale that establishes the 1/7(0;, 
scaling is still given by te independently of the initial inter¬ 
action strength; it is true for even an initially weakly inter¬ 
acting gas with 7(0;, 0) <C 1 and can emerge long before the 
instantaneous value of the interaction strength itself becomes 
“fermionic”, 7(x, t) ^ 1, due to its own scaling of j{x, t)(xt. 

In summary, we have analyzed the far-from-equilibrium dy¬ 
namics of the Lieb-Liniger gas in a quantum explosion sce¬ 
nario of a sudden expansion from the confining trap poten¬ 
tial. Considering a general class of power-law traps, we have 
found the asymptotic density profiles and the momentum dis¬ 
tributions of the expanding clouds using the stationary phase 
and local density approximations. The expansion is generally 
not self-similar, except for the strongly and weakly interacting 
gases released from a quadratic trap for which our results are 
in agreement with the known hydrodynamic scaling solutions. 
In all cases, the expanding clouds lose their phase coherence 
and display fermionic density fluctuations (not accounted for 
by the hydrodynamic theory) on a time scale R R/cq by 
which the expansion becomes ballistic. The zero-temperature 
results presented here are qualitatively valid for ksT po, 
however, our approach can be easily generalized to a nonzero- 
temperature initial state using the Yang-Yang approach EH, 
in which case the role of the initial local quasimomentum 
distribution f{k,x;t = 0) will be taken by its temperature- 
dependent counterpart to be found as in Refs. Il25l l3^ from 
the solutions to the Yang-Yang integral equations. 

The authors acknowledge stimulating discussions with M. 
Pustilnik and I. Bouchoule, and support by the ARC Discov¬ 
ery Project DP140101763. 
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I. THE MANY-BODY WAVE FUNCTION 


The Hamiltonian time evolution of the A^-particle wave 
function '^{xi,X 2 , • ■ •, xn] t) describing the dynamics of the 
Lieb-Liniger gas in free space, after its sudden release at time 
t = 0 from the confining trap potential, can be written down 
using an expansion in eigenfunctions j (xi, 0:2, • ■ • > xn) 
of the uniform Lieb-Liniger model. 


T'(xi, X2 ,..., xn; t) = J dki... dkN 

X h{ki,. .. ,kN)^{kj} (xi,X2,... ,X7v) (SI) 


Here, the eigenfunctions ■ ■ ■ ,xn) are character¬ 

ized by N different quasimomenta {kj} = ki,k 2 , ■ ■ ■ jk^ 
with eigenenergies E{kj} = /2m, whereas 

the expansion coefficients 5(fci,..., (cat) depend on 
the initial trapping potential and are normalized to 
J dki.. .dkN\b{ki,... = L The eigenfunc¬ 
tions }(xi,..., xat) in the fundamental region 

Xi < X2 < ... < Xn can be written down in explicit 
form as a sum of N\ partial waves [1], 


corresponds to N particles entering the collision, with the 
faster particles being behind the slower particles in consec¬ 
utive order from left to right. The term corresponding to the 
permutation obtained from the identity by swapping the in¬ 
dices i and j describes the scattering of particles with quasi¬ 
momenta ki and kj and therefore acquires a scattering phase 
9{ki — kj) = 2%sm~^\b?{ki — kj)/mg\. Proceeding in a simi¬ 
lar fashion one generates all other terms in the expansion (S2); 
in particular, the maximally crossed permutation term with 
Pj = N—j+1 corresponds to N particles leaving the collision, 
with the faster particles having overtaken the slower ones. 

Rearranging the quasimomenta {/cq} = kq ^, ,..., kq^ 

does not change the physical situation; however, the wave- 
function (S2) acquires the factor (—1)'3. Therefore, in order to 
avoid double counting of the physical states in the A^-particle 
bosonic wavepacket (SI), while at the same time ensuring 
that it remains symmetric, one requires that the coefficients 
h{ki, ..., A:at) are antisymmetric functions of quasimomenta. 

Substituting Eq. (S2) into Eq. (SI) and using the antisym¬ 
metry property of the expansion coefficients b{ki ,..., fcAr) 
gives a particularly simple expression for the initial state of 
the trapped Lieb-Liniger gas. 


T'{fc^.}(xi,... ,XAr) = 


E(-i) 


P ie(kp^ , 


y/{27T)^Nl 


(S2) 


enumerated by permutations P. Each of the Nl terms in the 
sum is multiplied by the sign of the permutation defined to 
be (—1)^ = ±1 for an even/odd permutation, Le., a permu¬ 
tation which can be obtained from the identity permutation 
{kj} = (ci, ^ 2 ) • • ■) AiAf by an even/odd number of transposi¬ 
tions of adjacent elements. The phase 


0(fci 




1 -I- i 


td(ki-k,) 


mg 


— I 


h?(ki-kj) 


= ^ tan ' 
3<l 


j<l - mg 

^b?{ki - kj)' 


mg 


(S3) 


arises from two-body scattering by the diffractionless delta- 
function potential of strength g (see main text). 

The structure of the eigenstates in Eq. (S2) has the follow¬ 
ing physical interpretation; in the fundamental region Xi < 
X2 < ... <XAr, choosing the decreasing arrangement of quasi¬ 
momenta fci > A:2 > . • ■ > Ajat in Eq. (S2) implies that the 
identity permutation term, in which Pj = j (j = 1,2,... N), 


T'(xi, X2 ,..., Xjv; 0) = f^ 2 TT)N /2 J dki...dkN 

xb{ki, ...,kN) (S4) 

valid in the fundamental region Xi < X2 < ... < xat. In any 
other region the symmetry of the bosonic functions under per¬ 
mutation of particles should be used. Equation (S4) describes 
the initial AV-body wavepacket; its time evolution after a sud¬ 
den removal of the confining potential is governed by the free 
Lieb-Liniger Hamiltonian and is given by 

T'(xi,X 2, ...,XN;t)= j dki...dkN 

Xb{ki, ...,kN) ^55^ 

which is the same as Eq. (1) of the main text. 

To evaluate the integral in Eq. (S5) using the stationary 
phase approximation, we note that the main contribution to 
the integral comes from the stationary phase points k* = 
mxj/ht corresponding to d(/>{kj)/dkj = 0, where (/{kj) = 
kxj — hkjt/2m. The integral is then calculated by approxi¬ 
mating the slowly varying function b{ki ,..., 
by its value at the stationary points and pulling it outside 
the integral; the remaining integral over the exponentials of 
phases (/{kj) are evaluated by using the Taylor expansion near 
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the points k* and keeping only the first two nonzero terms, 
(j){kj)~(f)[k*) + ^(j)"[k*){kj — k*Y. As a result, we obtain 


^00(^15^25 ■ ■ ■ 1 ^ 


/mx^/2 

rmxi 

mX]x[\ 

\nt) ' 

^ ht ’ 

■’ ht ) 


. ^ ) +* w j g-*7rtV/4 ^ 


which coincides with Eq. (4) of the main text. 


II. TWO-BODY CORRELATION FUNCTION 


interactions. Such scaling solutions can be obtained from the 
hydrodynamic approach [2-5], which we would like to illus¬ 
trate here for instruction. 

In the hydrodynamic approach, the evolution of the density 
and velocity fields, p{x, t) and u(x, t), under the variations of 
the trap frequency a;(f), is governed by the following coupled 
equations (see, e.g., [4, 5]), 

dtp + dx{pv) = 0 , 

dtv + vdxV = -dx ( )-(S14) 

\2 / mp 


The calculation of the two-body local (same point) correla¬ 
tion function. 


X J dxz . ■ -dXN \'i’oc{x,X,X3, ... ,xn; t)f , 


(SIO) 


requires the knowledge of the wavefunction at two coinciding 
points. The asymptotic wavefunction Eq. (S6) must, in this 
case, be modified to take into account the zero of the prefactor 

G(fci, fc2, .■.,kN) = b{ki,k2, ... (Sll) 


at ki = k 2 = mx/ht. Expanding this prefactor around the 
stationary points k* = mxj/ht to the second order, one ob¬ 
tains the required asymptotic wavefunction. 


'i>oo{x,X,X3, . . .,XN;t) = 
N 

i=i 


m 

iht \Ht) 
d^G\ 



^ — ITTNjA 


(S12) 


This form of the wavefunction manifests the effective 
‘fermionization’ via the additional factor 1 jt as compared to 
Eq. (S6). Substituting Eq. (S12) into Eq. (SIO) and using a 
variable change = mXj/Ht for j = 3 , ... ,N results in the 
following scaling of the pair correlation function with time: 
Gk^\x,x]t) oc 1/f^. Normalizing this with the product of 
densities pao{x,t), yields 


9^^\x,t) 


G^‘^\x,x;t) 1 

Poo{x,tY 


(SI 3) 


which is the result discussed in the main text, indicating the 
‘dynamical fermionization’ of expanding bosons. 


III. COMPARISON WITH THE HYDRODYNAMIC 
APPROACH 


As shown in the main text, the free expansion of the Lieb- 
Liniger gas from a quadratic {v = 2) trap leads to self-similar 
scaling solutions in the limits of very weak and very strong 


where P is the pressure. 

In the weakly interacting regime (70 <C 1 ), the pressure is 
given by P{x^ t) = gp{x, which follows from the equa¬ 
tion of state p, = gp and the Gibbs-Duhem thermodynamic 
relation p = [dP/dp)T. Using this and substituting the scal¬ 
ing ansatz for p{x, t) from Eq. ( 15 ) of the main text for ly = 2 
into the above hydrodynamic equations, yields, for uj{t) = 0 
{i.e., for complete opening of the confining trap), the follow¬ 
ing equation for the scaling parameter A(f): 

\ = ujI/X^. (S 15 ) 

The initial conditions are A( 0 ) = 1 and A( 0 ) = 0 . Mul¬ 
tiplying Eq. (S 15 ) by A one obtains the integral of motion 
A^ + 21x3 jX = const = 2ujq, which gives the asymptotic so¬ 
lution X{t) = \/ 2 (jjQt for t —)■ 00, i.e., the same result as the 
one discussed in the main text. 

In the Tonks-Girardeau regime (70 —> 00), the pressure is 
given by P{x,t) = p{x,t)^ /3m, from the equation of 

state p = p'^ /2m. In this case, substituting the scaling 

ansatz ( 19 ) of the main text (for v = 2) into the hydrodynamic 
equations with Lu(t) = 0 gives the following equation for the 
scaling parameter: 

X = ujI/X^. (S 16 ) 

The integral of motion is now A^ + Wq/A^ = const = Wg, 
which gives X{t) = wgf for t —> 00, in agreement with what 
we obtained in the main text. This result can also be obtained 
[6] within the exact treatment of the Tonks-Girardeau gas 
using time-dependent Bose-Eermi mapping and the scaling 
transformation applicable to the time-dependent Schrodinger 
equation for the single-particle orbitals. 
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